setwd("~/Dropbox/Gender Quota in Korea/_Pink-CollarWorker/Data")
library(foreign)
library(readstata13)
library(tidyverse)
library(cowplot)
library(cjoint);library(FindIt)

# Figure 1 #
# (1) Competency ####
dem.resp <- read.qualtrics("Working-Class+Women+Conjoint_October+13,+2021_13.12.csv",
                           responses = c("Q5.7", "Q6.7", "Q7.7"),
                           covariates = c("Q11.6", "Q11.8", "Pink.Collar", "Working.Class"),
                           respondentID = "ResponseId", new.format = TRUE)

rep.resp <- read.qualtrics("Working-Class+Women+Conjoint_October+13,+2021_13.12.csv",
                           responses = c("Q8.7", "Q9.7", "Q10.7"),
                           covariates = c("Q11.6", "Q11.8", "Pink.Collar", "Working.Class"),
                           respondentID = "ResponseId", new.format = TRUE)

data <- rbind(dem.resp, rep.resp)
data <- data %>% drop_na(selected)

data <- data %>% mutate(Education = factor(Education, levels = c("High school graduate", "College graduate", "Post-graduate (Master or Doctoral degree)")),
                        # Education = fct_relevel(Education, "High school graduate"),
                        Gender = fct_relevel(Gender, "Male"),
                        Working = factor(case_when(Former.Occupation == "Janitor" ~ "Working Class",
                                                   Former.Occupation == "Retail clerk" ~ "Working Class",
                                                   Former.Occupation == "Lawyer" ~ "White Collar",
                                                   Former.Occupation == "Business entrepreneur" ~ "White Collar")),
                        Working = fct_relevel(Working, "White Collar"))
model <- amce(selected ~  Age + Education + Working + Gender + Years.in.Politics, data = data, cluster=TRUE, respondent.id = "Response.ID")

labels <- c("Gender", "    Male", "    Female",
            "Occupational Background", "    White Collar","    Working Class",
            "Age", "   30", "   40", "   50", "   60", 
            "Education", "    High school graduate", "    University graduate", "    Post-graduate (Master or Doctoral degree)", 
            "Years in Politics", "    0", "    1", "    3", "    8"
)


estimates <- c(NA, 0, model$estimates[3][[1]][[1]], # gender
               NA, 0, model$estimates[4][[1]][[1]], # class
               NA, 0, model$estimates[1][[1]][[1]], model$estimates[1][[1]][[3]],model$estimates[1][[1]][[5]],
               NA, 0, model$estimates[2][[1]][[1]], model$estimates[2][[1]][[3]],
               NA, 0, model$estimates[5][[1]][[1]],  model$estimates[5][[1]][[3]],  model$estimates[5][[1]][[5]]
)

se <- c(NA, 0, model$estimates[3][[1]][[2]],
        NA, 0, model$estimates[4][[1]][[2]],
        NA, 0, model$estimates[1][[1]][[2]], model$estimates[1][[1]][[4]],model$estimates[1][[1]][[6]],
        NA, 0, model$estimates[2][[1]][[2]], model$estimates[2][[1]][[4]],
        NA, 0, model$estimates[5][[1]][[2]],  model$estimates[5][[1]][[4]],  model$estimates[5][[1]][[6]]
)

face <- c("bold", "plain", "plain",
          "bold", "plain", "plain",
          "bold", "plain", "plain", "plain", "plain",
          "bold", "plain", "plain", "plain",
          "bold", "plain", "plain", "plain", "plain")

type = c("gender", "gender", "gender",
         "occupation", "occupation", "occupation", 
         "age", "age","age","age", "age",
         "education", "education", "education", "education",
         "year", "year", "year", "year", "year")

issue <- data.frame(labels = factor(labels, levels = rev(labels)), 
                    estimates = estimates,
                    se = se,
                    face = face,
                    type = type)

issue <- issue %>% mutate(ul = estimates + 1.96 * se,
                          ll = estimates - 1.96 * se,
                          ref = ifelse(estimates == 0 | is.na(estimates), "reference", "normal"),
                          dv = "Competency")

issue
# (2) Candidate Support ####
dem.resp <- read.qualtrics("Working-Class+Women+Conjoint_October+13,+2021_13.12.csv",
                           responses = c("Q5.3", "Q6.3", "Q7.3"),
                           #                           covariates = c("Q11.6", "Q11.8", "Pink.Collar", "Working.Class"),
                           respondentID = "ResponseId", new.format = TRUE)
dem.resp <- dem.resp %>% drop_na(selected)
rep.resp <- read.qualtrics("Working-Class+Women+Conjoint_October+13,+2021_13.12.csv",
                           responses = c("Q8.3", "Q9.3", "Q10.3"),
                           #                          covariates = c("Q11.6", "Q11.8", "Pink.Collar", "Working.Class"),
                           respondentID = "ResponseId", new.format = TRUE)
rep.resp <- rep.resp %>% drop_na(selected)


data <- rbind(dem.resp, rep.resp)

data <- data %>% drop_na(selected)

data <- data %>% mutate(Education = factor(Education, levels = c("High school graduate", "College graduate", "Post-graduate (Master or Doctoral degree)")),
                        Gender = fct_relevel(Gender, "Male"),
                        Working = factor(case_when(Former.Occupation == "Janitor" ~ "Working Class",
                                                   Former.Occupation == "Retail clerk" ~ "Working Class",
                                                   Former.Occupation == "Lawyer" ~ "White Collar",
                                                   Former.Occupation == "Business entrepreneur" ~ "White Collar")),
                        Working = fct_relevel(Working, "White Collar"))
model <- amce(selected ~  Age + Education + Working + Gender + Years.in.Politics, data = data, cluster=TRUE, respondent.id = "Response.ID")
model$estimates

labels <- c("Gender", "    Male", "    Female",
            "Occupational Background", "    White Collar","    Working Class",
            "Age", "   30", "   40", "   50", "   60", 
            "Education", "    High school graduate", "    University graduate", "    Post-graduate (Master or Doctoral degree)", 
            "Years in Politics", "    0", "    1", "    3", "    8"
)


estimates <- c(NA, 0, model$estimates[3][[1]][[1]], # gender
               NA, 0, model$estimates[4][[1]][[1]], # class
               NA, 0, model$estimates[1][[1]][[1]], model$estimates[1][[1]][[3]],model$estimates[1][[1]][[5]],
               NA, 0, model$estimates[2][[1]][[1]], model$estimates[2][[1]][[3]],
               NA, 0, model$estimates[5][[1]][[1]],  model$estimates[5][[1]][[3]],  model$estimates[5][[1]][[5]]
)

se <- c(NA, 0, model$estimates[3][[1]][[2]],
        NA, 0, model$estimates[4][[1]][[2]],
        NA, 0, model$estimates[1][[1]][[2]], model$estimates[1][[1]][[4]],model$estimates[1][[1]][[6]],
        NA, 0, model$estimates[2][[1]][[2]], model$estimates[2][[1]][[4]],
        NA, 0, model$estimates[5][[1]][[2]],  model$estimates[5][[1]][[4]],  model$estimates[5][[1]][[6]]
)

face <- c("bold", "plain", "plain",
          "bold", "plain", "plain",
          "bold", "plain", "plain", "plain", "plain",
          "bold", "plain", "plain", "plain",
          "bold", "plain", "plain", "plain", "plain")

type = c("gender", "gender", "gender",
         "occupation", "occupation", "occupation", 
         "age", "age","age","age", "age",
         "education", "education", "education", "education",
         "year", "year", "year", "year", "year")

support <- data.frame(labels = factor(labels, levels = rev(labels)), 
                      estimates = estimates,
                      se = se,
                      face = face,
                      type = type)

support <- support %>% mutate(ul = estimates + 1.96 * se,
                              ll = estimates - 1.96 * se,
                              ref = ifelse(estimates == 0 | is.na(estimates), "reference", "normal"),
                              dv = "Candidate Preference")

support

results <- issue %>%
  bind_rows(support)
results

results <- results %>%
  mutate(dv = fct_relevel(dv, "Competency", "Candidate Preference"))
results

# AMCE plot:
ggplot(results, aes(x=labels, y=estimates, ymin=ll, ymax=ul, shape = type)) +
  labs(y = "Effect on Pr(Being Selected)") +
  geom_pointrange() + 
  theme_bw()  +
  #  coord_cartesian() + 
  coord_flip(ylim  = c(-0.3, 0.3)) + 
  scale_x_discrete(name="",labels=rev(labels)) +
  facet_wrap(vars(dv), nrow = 1) +
  geom_hline(yintercept = 0,colour="darkgrey",linetype=1) +
  #theme(axis.text.y = element_text(hjust=0)) +
  theme(
    axis.text.y = element_text(colour = "black",
                               hjust = 0 , vjust=.5,
                               face = rev(face), size = 14),
    axis.ticks.y = element_blank(),
    strip.text = element_text(size = 14),
    #axis.ticks.x = element_blank(),
    # axis.title.y = element_text(size = 7,angle=90,
    #                             vjust=.01,hjust=.1),
    legend.position = "none"
  ) + scale_shape_manual(values = c(18, 2, 20, 15, 7 )) 



# Figure 2: AMIE #
# (1) Competency ####
dem.resp <- read.qualtrics("Working-Class+Women+Conjoint_October+13,+2021_13.12.csv",
                           responses = c("Q5.7", "Q6.7", "Q7.7"),
                           covariates = c("Q11.6", "Q11.8", "Pink.Collar", "Working.Class"),
                           respondentID = "ResponseId", new.format = TRUE)

rep.resp <- read.qualtrics("Working-Class+Women+Conjoint_October+13,+2021_13.12.csv",
                           responses = c("Q8.7", "Q9.7", "Q10.7"),
                           covariates = c("Q11.6", "Q11.8", "Pink.Collar", "Working.Class"),
                           respondentID = "ResponseId", new.format = TRUE)

data <- rbind(dem.resp, rep.resp)
data <- data %>% drop_na(selected)
data <- data %>% mutate(Education = factor(Education, levels = c("High school graduate", "College graduate", "Post-graduate (Master or Doctoral degree)")),
                        Gender = fct_relevel(Gender, "Male"),
                        Working = factor(case_when(Former.Occupation == "Janitor" ~ "Working Class",
                                                   Former.Occupation == "Retail clerk" ~ "Working Class",
                                                   Former.Occupation == "Lawyer" ~ "White Collar",
                                                   Former.Occupation == "Business entrepreneur" ~ "White Collar")),
                        Working = fct_relevel(Working, "White Collar"))

data <- data %>% mutate(pairid = paste(respondent, task, sep="_"))

fit1 <- CausalANOVA(formula = selected ~  Age + Education + Working + Years.in.Politics + Gender , 
                    int2.formula = ~Working:Gender, data=data, pair.id = data$pairid, diff = TRUE,
                    cluster = data$Response.ID, nway=2)

cond.estimates<- ConditionalEffect(fit1, treat.fac="Working", cond.fac="Gender")$ConditionalEffects

labels <- c( "White Collar", "Working Class",
             "White Collar", "Working Class")



estimates <- c(cond.estimates[[1]][1,1],cond.estimates[[1]][2,1], 
               cond.estimates[[2]][1,1],cond.estimates[[2]][2,1])

ll <- c(cond.estimates[[1]][1,3],cond.estimates[[1]][2,3], 
        cond.estimates[[2]][1,3],cond.estimates[[2]][2,3])

ul <- c(cond.estimates[[1]][1,4],cond.estimates[[1]][2,4], 
        cond.estimates[[2]][1,4],cond.estimates[[2]][2,4])

condition <- c("Candidate Gender = Male", "Candidate Gender = Male", 
               "Candidate Gender = Female", "Candidate Gender = Female")
type = c("male", "male", 
         "female", "female")

issue <- data.frame(labels = factor(labels), 
                    estimates = estimates,
                    ll=ll,
                    ul=ul,
                    type=type, 
                    condition=condition)

issue <- issue %>% mutate(ref = ifelse(estimates == 0, "reference", "normal"),
                          dv = "Issue")


p1 <- ggplot(issue, aes(x=labels, y=estimates, ymin=ll, ymax=ul, shape=ref, color = type)) +
  scale_shape_manual(values = c(19,1)) +
  labs(y = "Effect on Pr(Being Selected)") +
  geom_pointrange() + 
  theme_bw()  +
  coord_flip(ylim = c(-0.2, 0.2)) + 
  scale_x_discrete(name="",labels=labels) + ylim(c(-0.15, 0.1)) + 
  facet_wrap(vars(condition), nrow = 2) +
  geom_hline(yintercept = 0,colour="darkgrey",linetype=1) +
  #theme(axis.text.y = element_text(hjust=0)) +
  theme(
    axis.text.y = element_text(colour = "black",
                               hjust = 0 , vjust=.5,
                               size = 14),
    axis.ticks.y = element_blank(),
    strip.text = element_text(size = 14),
    #axis.ticks.x = element_blank(),
    # axis.title.y = element_text(size = 7,angle=90,
    #                             vjust=.01,hjust=.1),
    legend.position = "none"
  )  + scale_color_manual(values = c("#e41a1c","#377eb8")) + ggtitle("Competency")
p1


# (2) Candidate Support ####
dem.resp <- read.qualtrics("Working-Class+Women+Conjoint_October+13,+2021_13.12.csv",
                           responses = c("Q5.3", "Q6.3", "Q7.3"),
                           covariates = c("Q11.6", "Q11.8", "Pink.Collar", "Working.Class"),
                           respondentID = "ResponseId", new.format = TRUE)

rep.resp <- read.qualtrics("Working-Class+Women+Conjoint_October+13,+2021_13.12.csv",
                           responses = c("Q8.3", "Q9.3", "Q10.3"),
                           covariates = c("Q11.6", "Q11.8", "Pink.Collar", "Working.Class"),
                           respondentID = "ResponseId", new.format = TRUE)

data <- rbind(dem.resp, rep.resp)

data <- data %>% drop_na(selected)
data <- data %>% mutate(Education = factor(Education, levels = c("High school graduate", "College graduate", "Post-graduate (Master or Doctoral degree)")),
                        Gender = fct_relevel(Gender, "Male"),
                        Working = factor(case_when(Former.Occupation == "Janitor" ~ "Working Class",
                                                   Former.Occupation == "Retail clerk" ~ "Working Class",
                                                   Former.Occupation == "Lawyer" ~ "White Collar",
                                                   Former.Occupation == "Business entrepreneur" ~ "White Collar")),
                        Working = fct_relevel(Working, "White Collar"))

data <- data %>% mutate(pairid = paste(respondent, task, sep="_"))


fit1 <- CausalANOVA(formula = selected ~  Age + Education + Working + Years.in.Politics + Gender , 
                    int2.formula = ~Working:Gender, data=data, pair.id = data$pairid, diff = TRUE,
                    cluster = data$Response.ID, nway=2)


ConditionalEffect(fit1, treat.fac="Working", cond.fac="Gender")

labels <- c( "White Collar", "Working Class",
             "White Collar", "Working Class")

cond.estimates<- ConditionalEffect(fit1, treat.fac="Working", cond.fac="Gender")$ConditionalEffects
estimates <- c(cond.estimates[[1]][1,1],cond.estimates[[1]][2,1], 
               cond.estimates[[2]][1,1],cond.estimates[[2]][2,1])

ll <- c(cond.estimates[[1]][1,3],cond.estimates[[1]][2,3], 
        cond.estimates[[2]][1,3],cond.estimates[[2]][2,3])

ul <- c(cond.estimates[[1]][1,4],cond.estimates[[1]][2,4], 
        cond.estimates[[2]][1,4],cond.estimates[[2]][2,4])


condition <- c("Candidate Gender = Male", "Candidate Gender = Male", 
               "Candidate Gender = Female", "Candidate Gender = Female")
type = c("male", "male", 
         "female", "female")

support <- data.frame(labels = factor(labels), 
                      estimates = estimates,
                      ll=ll,
                      ul=ul,
                      type=type, 
                      condition=condition)

support <- support %>% mutate(ref = ifelse(estimates == 0, "reference", "normal"),
                              dv = "Candidate Preference")


p2 <- ggplot(support, aes(x=labels, y=estimates, ymin=ll, ymax=ul, shape=ref, color = type)) +
  scale_shape_manual(values = c(19,1)) +
  labs(y = "Effect on Pr(Being Selected)") + 
  coord_flip(ylim = c(-0.2, 0.2)) + 
  geom_pointrange() + 
  theme_bw()  +
  scale_x_discrete(name="",labels=labels) + ylim(c(-0.15, 0.1)) + 
  facet_wrap(vars(condition), nrow = 2) +
  geom_hline(yintercept = 0,colour="darkgrey",linetype=1) +
  #theme(axis.text.y = element_text(hjust=0)) +
  theme(
    axis.text.y = element_text(colour = "black",
                               hjust = 0 , vjust=.5,
                               size = 14),
    axis.ticks.y = element_blank(),
    strip.text = element_text(size = 14),
    #axis.ticks.x = element_blank(),
    # axis.title.y = element_text(size = 7,angle=90,
    #                             vjust=.01,hjust=.1),
    legend.position = "none"
  )  + scale_color_manual(values = c("#e41a1c","#377eb8")) + ggtitle("Candidate Preference")
p2


library(ggpubr)
ggarrange(p1, p2, ncol=2)

